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ABSTRACT 

We simulate the formation of a low metallicity (10“^ Zq) stellar cluster in a dwarf galaxy at 
redshift ^ ^ 14. Beginning with cosmological initial conditions, the simulation utilizes adap¬ 
tive mesh refinement and sink particles to follow the collapse and evolution of gas past the 
opacity limit for fragmentation, thus resolving the formation of individual protostellar cores. 
A time- and location-dependent protostellar radiation field, which heats the gas by absorption 
on dust, is computed by integration of protostellar evolutionary tracks with the MESA code. 
The simulation also includes a robust non-equilibrium chemical network that self-consistently 
treats gas thermodynamics and dust-gas coupling. The system is evolved for 18 kyr after the 
first protostellar source has formed. In this time span, 30 sink particles representing protostel¬ 
lar cores form with a total mass of 81 Mq. Their masses range from ^0.1 Mq to 14.4 Mq 
with a median mass ^0.5 — 1 Mq. Massive protostars grow by competitive accretion while 
lower-mass protostars are stunted in growth by close encounters and many-body ejections. In 
the regime explored here, the characteristic mass scale is determined by the temperature fioor 
set by the cosmic microwave background and by the onset of efficient dust-gas coupling. It 
seems unlikely that host galaxies of the first bursts of metal-enriched star formation will be 
detectable with the James Webb Space Telescope or other next-generation infrared observa¬ 
tories. Instead, the most promising access route to the dawn of cosmic star formation may lie 
in the scrutiny of metal-poor, ancient stellar populations in the Galactic neighborhood. The 
observable targets that correspond to the system simulated here are ultra-faint dwarf satellite 
galaxies such as Bootes II, Segue I and II, and Willman I. 
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1 INTRODUCTION 

The cosmic ‘dark ages’ was a transformative epoch in the history 
of the Universe, witnessing the birth of the first stars, galaxies, and 
supermassive black holes, the beginning of reionization, and the 
onset of chemical enrichment ( |Barkana & Loeb|2001| ). Despite its 
importance, it is also one of the least understood epochs. Metal- 
poor stars and dwarf galaxies, likely relics from the high-redshift 
Universe, present an opportunity to use local observations to inves¬ 
tigate the earliest stellar generations and their role in driving early 


cosmic evolution (e.g.,|Beers & Christlieb|2005| 

Frebel & Bromm| 

|2012||Karlsson et al.|2013||Frebel & Norris|2013 

1 . To capitalize on 


the power of stellar and the Galactic dwarf satellite galaxy ‘archae¬ 
ology’ as a probe of primeval stellar populations, it is necessary 
to develop a basic theoretical understanding of star formation as 
it occurred in the first galaxies at high cosmic redshifts and low 
metallicites. 

Modeling star formation in the first galaxies is challenging in 


part because we still lack a robust, predictive theory of present-day 


star formation (see reviews by|McKee & Ostrikerj 

2007 

,|Zinnecker| 

|& Yorke||2007[|Krumholz et al.||2014[|Tan et all 

|2014 

1 , let alone 


for that occurring in the early Universe ( |Bromm|2013| l. Neverthe¬ 
less, one plausible conjecture is that gas metallicity, likely highly 


sub-solar in the first galaxies (e.g., [Wise & Abelp008||Greif et^ 
2010|), moder ated the thermodynamics of star forming gas (e.g., 
Larson] 1985| ). 

The thermodynamic evolution of collapsing gas concentra¬ 
tions, determined through the competition of heating and cool¬ 
ing processes, is known to be crucial for modulating fragmenta¬ 
tion (e.g., |Li et al. |200^ and for potentially fixing the character¬ 
istic (i.e., median) mass of the stellar initial mass function (IMF; 
e.g., jMasunaga et ^|1998[ |0mukai||2000[ |Larson||2005[ jJappsenj 

jet al.|2005| l. In preset-day molecular clouds, the primary heating 

sources are dust grain photoelectric absorption and Coulomb colli¬ 
sions with cosmic rays. Cooling occurs through line emission (C^ 
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or CO) and dust grain thermal emission, at low and high densi¬ 
ties, respectively. While the precise balance of these processes de¬ 
pends largely on the metallicity, column density, and strength of 
external radiation fields, star forming gas evolves approximately 
isothermally until the column density has increased enough for the 
dust to become opaque to its own cooling radiation. After that, the 
evolution is quasi-adiabatic. 

In Population III (Pop III) stars which form from metal-free 
gas, the thermodynamic evolution is markedly different, dictated 
instead by the physics of molecular hydrogen, H 2 , the most effec¬ 
tive coolant in low-temperature (T < 10^ K), metal-fee gas (e.g., 
SiIklfT^ [Abel et al.|[2000l |Bromm et al.1[200^ [Yoshida et alT] 


If metal enriched star formation promptly followed Pop III star 
formation ( e.g., [Ricotti et al.|2002| [2008[[Ritter et al.|20T^|Whale'n| 


et al. 


2013t , Mvir 10^ Mq ‘minihalos’ would have been the host 


2006| l. This is believed to make the characteristic Pop III stellar 
larger (e.g., Stacy et al.|20I0[|Clark et al.|20II[[Greif et al.| 


|20II[ |Hirano et al.|20I4|l than t he typical Galactic stellar mass of 


0.1 Mq (e.g.,|Chabrier|2003ll. 


The salient heating and cooling processes at temperatures be¬ 
low ~ 10^ K are sensitive to metallicity. Using radiation hydrody- 
namical (RHD) simulations, [Myers et al.| ( [20TT] ) and |Bate| ( |2014| ) 
examined the effect of metallicity variation on the star formation 


sites of the very first Pop II stellar clusters. A minihalo, however, 
does not meet other intuitive criteria for a bona fide galaxy, such as 
the requirement that it be able to retain photoionized gas. This tax¬ 
onomic subtlety aside, the physical state of primordial star forming 
systems just prior to the onset of metal-enriched star formation has 
been explored in a number of studies, all reaching largely compat¬ 
ible conclusions. This physical state defines the initial conditions 
for metal-enriched star formation in the first galaxies. 

Simulations of the transformation of a gas cloud into stars are 
traditionally initialized with an isolated spherical or quasi-spherical 
overdensity. A spectrum of random velocity perturbations is im¬ 
posed at initialization, or hydrodynamic turbulence is externally 
driven until saturation. The parameters of the initial cloud and 
of the velocity (or perturbing force) fiuctuations are often chosen 
ad hoc, reflecting a focus on the rudimentary mechanisms of star 
formation, rather than on making quantitative predictions in var¬ 
ied, realistic, complex astrophysical settings. Properties of the re- 


process in the context of present-day star formation. They found 

suiting stellar systems, such as their IMFs {|Klessen|200I 

Martel 

that down to 10“^ Z©, metallicity had little effect on the IMF and 

et al. 

20061 [Clark et al.|2008| [Urban et al.|20I0l[Girichi( 

is et al. 

other properties of simulated stellar systems, in agreement with 

2011 

1 , or the mode of the stellar mass fixation (e.g., ‘core accretion’ 


though, dust and gas were assumed to be perfectly collisionally 
coupled, so that the sole effect of lowering metallicity was a pro¬ 
portional reduction in dust opacity. In reality, the assumption of 
perfect dust-gas coupling breaks at low densities and metallicities, 
an effect that can significantly alter the gas cooling rate and ther¬ 
modynamic behavior. 

The impact of metals on the thermodynamic evolution of col¬ 
lapsing gas concentrations has also been studied in the theoret¬ 
ical context of the Pop III to Pop II star formation mode tran¬ 
sition. Idealized, one-zone thermodynamic models generally sug¬ 
gest that once the gas metallicity has exceeded some critical value 
Zcrit, the thermodynamical behavior, and thus characteristic frag¬ 
mentation mass, change dr astically, either due to metal line cooling 

(■^crit 


10 “ 


^o;e.g., |Bromm & Loci 


fi]2003|p 

[2006[[Safranek-Shrader et al.|2010[|Omukai et al.|2010| l, or due to 

the onset of gas-dust collisional coupling (Zcrit ~ Zq \ e.g.. 


Santoro & Shull 


[nell et al.|2006| ), exhibit sensitivity to the model parameters. En¬ 
couragingly, state-of-the-art radiation-hydrodynamical (RHD) sim¬ 
ulations have successfully delivered stellar clusters with conver¬ 
gent IMFs and other statistical measures matching observed low- 
luminosity star forming regions such as the Orion nebula cluster 
(e.g., [Krumholz et al.|2012[[Bate|2012[ l, though fine tuning of cer¬ 
tain parameters is still required. Star formation simulations that be¬ 
gin from ab initio initial conditions—either seeded from global 
galactic disk simulations, or, ultimately, by infiationary random 
fields—are needed to settle these uncertainties. 

In the first simulation to deliver low-mass protostars directly 
from cosmological random fields { [Safranek-Shrader et al.|20I4^ , 
we simulated the assembly of a metal-poor (10“^ Zq) stellar sys- 
tem in a redshift 2 : ~ 14 atomic cooling halo. The initial conditions 
were excised from the cosmological box of [Safranek-Shrader et al.[ 
P014b[ ) evolved to peak densities in gravitationally collapsing gas 
Omukai et al.[[2005| [Schneider et al.[[2006[ [Schneider & Omukai[ clumps of ~ 10^ cm“^ . The excised box was centered on a single 

dense clump that formed via metal-line-cooling-induced thermal 
instability. We adaptively refined the computational grid and intro¬ 
duced Lagrangian sink particles to track gravitational collapse of 
gas to densities ~ 10^^ cm“^ and resolved the formation of in¬ 
dividual protostellar cores. Over the maximum simulation time of 
7 kyr permitted by computing resources, 40 protostellar cores 
formed with final masses ranging from Mq to 2.5 Mq. The 
stellar IMF above 0.1 Mq was tentatively shallower (top-heavier) 
than the Salpeter power law. Heating by protostellar accretion ra¬ 
diation did not have an impact on the low-mass IMF, but the time 
extent of the simulation following the onset of star formation, cor¬ 
responding to only a tenth of the local free fall time, was insufficient 
for massive, luminous protostars to form and appreciably heat the 
dust. In principle, we could have integrated the system as much as 
~ 10 times longer while still avoiding boundary artifacts from the 
box excision. Such an extension of the simulation would have al¬ 
lowed the protostellar mass function to grow to higher masses and 
much higher luminosities. 

In the present paper, we extend the simulation of [SafraneE^ 
[Shrader et al.[ ( [2014^ while crucially calibrating individual proto¬ 
stellar luminosities to evolutionary tracks computed with the new 


confirmed this ([Bromm et al. |200I| [Smith & Sigurdsson[ 

2007 

Clark et al.[|2008[ [Smith et al.| 

2009 1 [Dopcke et al.l 

2011 

2013 

Safranek-Shrader et al.|20I4a|b 

[Bovino et al.|20I4ll 

IMeece et al.[ 

2014|. 


A thorough understanding of star formation in the first galax¬ 
ies requires placing the process in its proper cosmological context. 
The first galaxies are anticipated to form in dark matter halos with 
virial temperatures Tvir ~ 10^ K c orresponding to viria l masses 
Mvir ~ 10® [(1 + «)/10]"®/^ Mq ^Oh & Haiman|2002||Bromm 
[& Yoshida|20II[ l. In these ‘atomic cooling halos’, energetic super¬ 
novae, infall of baryonic matter from the cosmic web, and perhaps 
even thermal instability, all contribute to the excitation of turbu- 


lence in the star-forming gas (ISM; [Wise & Abel|2007 

[Greif et al. 

2008 

[Wise et al.|20081 [Prieto et al.|20I21[Safranek-S] 

irader et al. 

2012 

to avt 

1 . Previous stellar generations could have polluted these halos 
jrage metallicities Z > Zcrit (e.g.,|Tomatore et al.[2007[[Wise| 

[& Abel|2008[[Greif et al.[20I0[[Maio et al.[20I0l[Ritter et al.[20I2[|, 


which would have enabled the formation of Pop II stars (but metal 
pollution io Z > Zcvit cannot be taken for granted, see [Ritter et al.[ 
|20T4] ). 
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Star Cluster Formation in the First Galaxies 3 


Stellar evolutionary package Modules for Experiments in Stellar As¬ 
trophysics (MESA; [Paxton et al.|2011|[2013| l. In addition to these 
improvements, we provide the technical detail left out in the Letter- 
length report of our early results in jSafranek-Shrader et aL| ( |2014a| t. 

We organize this paper as follows. In Section we describe 
the methodology and simulation setup. This includes a brief de¬ 
scription for the parent cosmological simulation from which our 
initial conditions were excised as well as a detailed description of 
our chemo-thermodynamic and proto stellar evolutionary model. In 
Section we describe the results of the simulation with a special 
emphasis on the formation and growth of sink particles and the 
thermodynamics of the gas and dust in the course of gravitational 
collapse. In Sectionj^we discuss the implications of our results and 
summarize our conclusions. 


2 METHODOLOGY 
2.1 Simulation setup 


We performed our simulation with the adaptive-mesh-refinement 
(AMR) hydrodynamics code FLASH ( jFryxell et al.|2000| l, version 
4. The initial conditions were extracted from the cosmological sim¬ 
ulation of jSafranek-Shrader et al.| { |?014b| ) and consisted of a cubi¬ 
cal region of size 0.52 pc containing a total gas mass of 390 Mq. 
In this parent simulation, a cosmological volume of 1 comoving 
Mpc^ was evolved until an atomic cooling halo with a virial tem¬ 
perature of Tvir ~ 10^ K formed in the box. Emulating a radiative 
background due to star formation outside the cosmological box, we 
introduced a global Lyman-Werner (LW) radiation field with inten¬ 
sity Jlw,2i = 100 which had the effect of photodissociating H 2 
in low-column-density gas and preventing Pop III star formation in 
the progenitor minihalos[^The atomic cooling halo virialized with 
mass Mvir ~ 2 x 10^ M© at z = 13.8. At that point, we set the gas 
metallicity inside its virial radius from zero to 10 modeling 

metal enrichment by preceding Pop III stars in a highly idealized 
fashion. The excision of a dense clump was performed when the 
peak density reached 10^ cm“^. The gravitational potential in the 
excised region was strongly baryon-dominated and so after exci¬ 
sion we neglected dark matter. The mass-weighted average density 
inside the box was 9 x 10^ cm“^ corresponding to a free-fall time 
of tff = (37r/32Gp)^/^ « 50kyr. 

The excised hydrodynamic data cube was used to refine the 
simulation to densities and time scales unattainable in the larger 
cosmological box while still retaining direct causal dependence 
on the cosmological random fluctuations. We imposed refiective 
boundary conditions crudely approximating a pressure-confined 
environment. To prevent contamination of our results by bound¬ 
ary effects, we could run the simulation only for much shorter than 
the sound crossing time from box edge to center, ~ 0.4 Myr. We 
always resolved the Jeans length 


Lj 



1/2 


( 1 ) 


by at least 24 grid cells. Through experimentation we found that 
resolving the Jeans length by fewer cells produced a significantly 


^ Here, Jlw,21 is the mean intensity at the center of the LW bands, 
12.4 eV, expressed in the units of 10“^^ ergs“^ cm“^ Hz“^ sr“^. This 
is closely related, thought not exactly equal, to J 21 , the radiation intensity 
at the Lyman limit. 


different turbulent morphology and non-convergent protostellar 
growth trends (see, also, |Federrath et al.|201 Ij and jTurk et al.|2012| |. 


2.2 Sink particles and the formation of individual stars 

Above some critical column density, collapsing gas concentrations 
generally transition from isothermal to adiabatic evolution when 
the continuum opacity exceeds unity and radiative cooling loses 
efficacy. This is the point when the energy input rate from gravi¬ 
tational compression exceeds the maximum allowed radiative loss 
rate. The thermal Jeans mass at the isothermal-to-adiabatic transi¬ 
tion is known as the opacity limit for fragmentation ( [Rees l|TW6l 
[Low & Lynden-Bell||1976| ). It represents the minimum mass of a 
gravitationally collapsing gas clump. jOimukaij ( [2000] ) showed that 
if the optical depth is estimated as the opacity multiplied with the 
local Jeans length, as is appropriate in self-gravitating collapsing 
gas, the density at which the optical depth equals unity is indepen¬ 
dent of opacity (and thus of metallicity) and is given by 

^ 12 T^CTgB^H ^ ^ 

where T is the gas temperature, mu is the hydrogen atom mass, and 
ku and ctsb are the Boltzmann and Stefan-Boltzmann constants, 
respectively. Thus, to resolve the formation of individual stars, we 
must follow the gas collapse past ricrit- 

Sink particles (hereafter sinks) are a commonly employed in 
simulations of star formation. Originally introduced by [Bate et aL] 
( |1995| ), they have been widely used in both Eulerian (Godunov- 
type) and Lagrangian (smoothed particle hydrodynamics; SPH) 
methods. As self-gravitating gas collapses and density increases, 
the Eulerian grid refines itself to resolve the local Jeans length by a 
minimum number of resolution elements. Eventually, the Courant- 
Eriedrichs-Lewy limited time step becomes prohibitively short. To 
continue the simulation past the initial runaway collapse, it be¬ 
comes necessary to introduce a subgrid model for collapsed sites. 
This is done by inserting Lagrangian, collisionless, gravitating par¬ 
ticles into sites where the density has exceeded a threshold. The 
particles accrete excess density gas from the grid in a momentum- 
conserving fashion, keeping the Jeans length in check and obviating 
further grid refinement. 

Sinks, however, are not just a computational shortcut. If 
used correctly, they can represent distinct, localized, gravitation¬ 
ally collapsed gas concentrations. They make it straightforward 
to record the mass accretion rate and other subgrid features of 
the histories of these concentrations. To regard sinks as individ¬ 
ual protostellar cores (and minimize the chance of unresolved sub¬ 
fragmentation occurring within a sink), we set the density at which 
sinks form to rigink = 4 x 10^^ cm“^, nearly an order of magni¬ 
tude higher than the estimate given by Equation g. Physically, this 
means that gas that reaches rigink becomes incorporated in a well- 
defined, pressure-supported, quasi-adiabatic core in which gravita¬ 
tional fragmentation is suppressed. 

In addition to the density threshold, for sink particle creation 
we also require that the fiow in the cell with n > rigink be con¬ 
verging V • V < 0, the gravitational potential be a local mini¬ 
mum, and a small control volume around the cell be gravitation¬ 
ally bound, E^therm + ^kin + E^grav < 0. Hcrc, E^kin is the total 
kinetic energy evaluated with respect to the center-of-mass veloc¬ 
ity. These additional checks are essential to prevent the formation 


^crit 
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of spurious sink particles that do not represent physical gravita¬ 
tionally collapsing sites ( |Federrath et al.|20T0] |. In each hydrody¬ 
namic time step, grid cells within the sink’s accretion radius of 
Ta.cc — 20 AU = 2.5 Axmin with densities of hydrogen nuclei 
riH > risink transfered a fraction (nn — nsink)/'^^H of their mass to 
the sink if the gas was gravitationally bound to the sink and had a 
velocity with a negative radial component relative to the sink. Here, 
Axmin is the cell size at the highest level of grid refinement. Sink 
particle motion is sub-cycled with a leapfrog integration scheme. 
Sinks are not allowed to merge with each other. For further details 
and limitations of the sink particle approach, see |Federrath et al.| 


2.3 Thermodynamical model 

Our thermodynamical model and non-equilibrium chemical net¬ 
work were as in [Safranek-Shrader et al.| ( |2010| |2012| |2014b| t, but 
now with the dust processes described in jOmukai et al.j ( |2005| ) 
and additional processes that become relevant at high densities, 
^ 10® cm“®. At sub-solar metallicities and low densities, dust 
and gas cannot be assumed to be thermally coupled. In general, 
their temperatures differ. We self-consistently calculated the dust 
temperature Td by assuming balance between cosmic microwave 
background (CMB) and protostellar radiation absorption, thermal 
emission, and inelastic gas-dust collisions. We numerically solved 
the balance relation 

4crsB(Td — TcMB)^d(Td)pPesc — -T- 

tcoll 

+ E(^) (3) 


for the dust temperature Td. Here, p is the gas density, Tcmb 
is the CMB temperature, ri is the distance to the zth sink, and 
Ttot,i is the total protostellar luminosity of the zth sink (see Sec- 
tion |2.4| ). We assumed that the number density of dust grains scaled 
linearly with metallicity rid = pD{ZjZ q)I rrid, where D is the 
Galactic dust-to-gas ratio, which we took to be T) = 0. 01, and 
rrid = 1.3 X 10“^^ g is the effective dust grain mass ( [Cazaux 


|& Spaan^|2004j ). The collision time between gas and dust par¬ 
ticles is — TinUdViif, where — 0.1 pm is the aver¬ 
age dust grain cross section ( jCazaux & Spaans||2004| l, vn is the 
average speed of hydrogen nuclei, and / 0.4 accounts for 

non-hydrogenic species ( [Schneider et al.|20'06| l. The Planck mean 


opacity of dust grains, /^d(Td), was taken from Semenov et al. 
( [2003 j ). We assumed that /^d(Td) scales linearly with metallicity 
and adopted a density-independent dust sublimation temperature 
of 1500 K. Thermal emission from dust grains is attenuated by a 
factor /3esc = min(l, appropriate in the regime of optically- 

thick radiative diffusion (e.g., [Masunaga et al.|19'^ . The contin¬ 
uum optical depth is given by Tcont = (^d + ^g)pTj, where Lj 
is used as a local estimate of the physical extent of a gravitation¬ 
ally collapsing core and /^g(p, T) is gas Planck mean opacity. We 
ignored the metal contribution to hig and took the gas opacity from 
[Mayer & Duschl[ ( [2005[ ). Ignoring metals in Kg is justified because 
the bulk of metal opacity is accounted for in Kd. 

To determine the metal fine structure line cooling rate, we uti¬ 
lized the same procedure as in [Safranek-Shrader et al.[ ( [2014b[ t, 
but now model radiation trapping by line and dust opacity. Since 
the level populations and line escape probabilities depend on each 
other, obtaining a self-consistent cooling rate requires iterative so¬ 
lution (e.g., [Takahashi et al.|1983[[0mukai|2000| ). We employed a 



Figure 1. Protostellar radius (top panel) and effective temperature (bottom) 
as a function of stellar mass from our MESA evolutionary tracks of accreting 
protostars. The colors refer to different values of the constant mass accretion 
rate (see legend). Deuterium burning is active at M^c ~ (0.3 — 0.5) Mq 
for M* = lO“®M 0 yr“^ and at M* ~ (0.5 — 2) Mq for M* = 
10“^ Mq yr“^ and is manifested in a shallow bump in the stellar ra¬ 
dius. Radius and the effective temperature show a notable, sudden increase 
with core contraction and the associated convective-to-radiative transition. 
Kelvin-Helmholtz contraction toward the zero-age main sequence (ZAMS) 
follows thereafter. Above ~ 20 Mq, the 10“® Mq yr“^ track begins as¬ 
cent of the red giant branch. 


local estimate of the Sobolev length, 

Tsob = 


|V-v| 


(4) 


to approximate the size of the shielding region. 

Our computations accounted for cooling by ro-vibrational 
lines of H 2 . Molecular hydrogen forms via the H“ intermediary 
and on the surfaces of dust grains. Above n ~ 10® cm“®, H 2 be¬ 
gins to form efficiently in three-body reactions. We track the la¬ 
tent heat of H 2 formation and dissociation. The reduction in the 
ro-vibrational H 2 cooling rate due to line trapping was modeled 
using an escape probability formulation based on the same local 
estimate. E quation 0- of the Sobolev length (see, e.g., [Yoshida[ 
et al.|2006 ( ^ Our model also included collisionally induced emis¬ 
sion (CIE) from H 2 , though this process does not become a signifi¬ 
cant cooling channel until densities exceed ~ 10^^ cm“®, a regime 
not explored here. 


2.4 Protostellar growth and feedback 

Radiation emitted from the surfaces and disks of accreting proto¬ 
stars can modify conditions in the ambient star-forming gas (e.g., 
[Offner et al.|2009[ t. The summand in the last term of Equation ^ 


^ This approach may overestimate the H 2 line emission escape fraction by 
a large factor in certain regimes jGreif[2014| . 
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Star Cluster Formation in the First Galaxies 5 


represents the heating rate of dust by the radiation of the zth proto¬ 
star located at distance with a total luminosity 


Ltot — L'lYit Lace 

that is a combination of the intrinsic luminosity 

Lint — 47ri?>^(7SBX'gff , 

and accretion luminosity 


= /a. 


GM.M. 


(5) 


( 6 ) 


(7) 


where face is a dimensionless coefficient quantifying radiative ef¬ 
ficiency of accretion. By treating sinks as sources of radiation, we 
identify with the sink mass and M* with rate of accretion onto 
the sink. The accretion rate was smoothed over a 2 yr window cor¬ 
responding to about 5 hydrodynamical time steps. We assume a 
fixed radiative efficiency of face = 0.75 (e.g., |Offner et al.|2009| ). 
This choice hides our ignorance of the detailed geometry of the 
protostellar accretion fiow, as we explain below. 

The radius of the protostellar photosphere and the effective 
temperature Tgff evolve significantly over the course of protostel¬ 
lar growth. We used MESA to tabulate the protostellar radius and 
effective temperature as a function of protostellar mass and accre¬ 
tion rate. Each MESA track was computed with a constant (time- 
independent) mass accretion rate . The entropy of the accreted 
gas was set to that of the stellar surface, approximating accretion 
via a protostellar disk. We assumed standard Big Bang nucleosyn¬ 
thesis (BBN) abundances for hydrogen and helium, a metallicity 
of 10 “^^ 0 , and the standard solar abundance pattern ( Grevesse 
|& Sauval||1998| ). For Tgff < 43,000 K, MESA photospheric ta¬ 
bles were used as the atmospheric boundary condition, while for 
Teff > 43, 000 K, an Eddington grey atmosphere was assumed. For 
deuterium burning we used the prescription of | Weiss et al.| ( |200'4] ) 
for convection with a mixing-length-to-pressure-scale-height ra¬ 
tio of qml = 2.0. The MESA integrations were initialized with 
M* = 0.1 Mq and R^ = 20i^o and stopped when either the 
stellar mass reached 50 Mq or the star left the main sequence. In 
Figure we show the protostellar radius and effective temperature 
as a function of protostellar mass for four different values of M*. 

To interpolate a sink particle’s R^ and Tgff from tabulated 
tracks, instead of the sink’s rapidly fluctuating instantaneous ac¬ 
cretion rate, we used its lifetime averaged accretion rate 


t) = 


M.jt) 

t tform 


( 8 ) 


where tform is the sink formation time. The interpolated radii and 
temperatures were then used in Equations ([^ and ([^ to compute 
the sink’s total luminosity. We neglected the luminosities of sinks 
with M* <0.1 Mq. This procedure produced stellar radii and ef¬ 
fective temperatures in good agreement with the low-metallicity, 
disk-accretion models of |Hosokawa et al.| ( |2010| ). 

The geometry of the accretion fiow can have a significant ef¬ 
fect on protostellar evolution (e.g., [Hartmann et al.|1997[[Behrend 
& Maederj[20011 [Hosokawa & Omukai[[2009[ [Hosokawa et al. 
2010| l. The model in which the accretion flow is spherically sym¬ 
metric and in free fall is known as ‘hot accretion’. The material ac¬ 
quires a high entropy upon passing an accretion shock at the stellar 
surface and blankets the surface. This produces bloated protostellar 
radii. If accretion is channeled through a circumstellar disk, the ma¬ 
terial loses entropy prior to joining the star and the blanketing effect 
is absent, so the protostellar radius is comparably smaller. Accre¬ 
tion collimated by the protostellar magnetosphere may belong in 


between these limits, receiving an entropy increment in the accre¬ 
tion shock, but with little blanketing. Thus protostellar evolution 
is not uniquely determined by the accretion rate but also depends 
on the accretion geometry. By simulating the evolution of accret¬ 
ing protostars on the Hertzsprung-Russell diagram, [Hosokawa et al.[ 
( [ 201 1 [ ) found that protostellar tracks with pure disk accretion sig¬ 
nificantly underpredict the luminosity of hot (T > 4000 K) proto¬ 
stars. 

Observationally, it is challenging to determine the accretion 
geometry due to the heavy obscuration of embedded protostellar 
sources. Disks are known to be common around massive stars (e.g., 
[Patel et al.|2005|[Cesaroni et al.|2007| ) and simulations of massive 
star formation typically show circularization of accretion streams 
into disks (e.g., [Yorke & Sonnhalter[[2002[ [Banerjee et al.[[200^ 

[Krumholz et al.|2007[[2009t . Thes^imulations, however, do not 

resolve the protostellar magnetosphere and cannot exclude the pos¬ 
sibility that the material accreting through a disk eventually chan¬ 
nels onto the magnetic poles. 

In [Safranek-Shrader et al.[ ( |^14a[ t we crudely modeled the ac¬ 
cretion luminosity and neglected the intrinsic luminosity. This sim¬ 
plification is valid for low-mass (< 3 Mq) protostars. Specifically, 
we estimated the accretion luminosity with Equation 0 after sub¬ 
stituting the [Stahler et al.[ ( p~986[ ) estimate of the protostellar radius 


— 2GRq 


( M. 

\Mq 


(9) 


The latter estimate, however, assumes spherically symmetric, hot 
accretion. In contrast, we computed the MESA tracks consistent 
with the premise that gas arrives at the protostellar surface through 
a cold, thin disk. Since hot accretion models lead to protostellar 
radii potentially orders of magnitude larger than cold accretion 
models (e.g., [Hosokawa et al.[[2010[ ), the sink particle accretion 
luminosity in this present study is much larger than in [Safranek-1 
[Shrader et al.[ ( [2014^ . 

Finally, we emphasize that our prescription for protostellar ra¬ 
diative feedback is entirely based on local approximations, mod¬ 
eling the effect of optical depth with the escape probability for¬ 
malism. Radiative transfer simulations like that of IKrumholz et al.1 
\2012\ are required to raise the modeling of primordial star forma¬ 
tion to a new level of realism. 


3 RESULTS 

3.1 Structure of the Collapsing Cloud 

The initial conditions, as extracted from the parent cosmological 
simulation, consist of an ellipsodial pre-stellar clump on the verge 
of gravitational collapse. The spherically averaged density profile 
is p oc from ^ 4 x 10^ AU to the edge of the computational 
box. At radii <10^ AU, the density profile is relatively shallower, 
and further inside 10^ AU, it levels off to nn ~ 10^ cm“^. Similar 
density profiles have been used by other groups to initialize isolated 
star-forming spherical clouds (e.g., [Krumholz et al.|2012[ > and are 
consistent with observations of massive star forming clumps (e.g., 
[Beuther et al.|2006] t. 

Figure 1^ shows snapshots of the morphological evolution of 
the central ~ 15, 000 AU from the initial ellipsoidal clump un¬ 
til 51 kyr later. The first sink forms 33kyr after the beginning 
of the simulation in gas compressed by a large-scale fiow conver¬ 
gence inherited from the parent cosmological simulation. In the 
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Figure 2. Mass-weighted line-of-sight projections of gas density in the x — y plane. The top-left panel shows the initial state of the gas as extracted from the 
parent cosmological simulation. Time increases from left to right and from top to bottom, with f = 0 shifted to the formation of the first sink. Circles mark 
projected sink locations with the circle size increasing with sink mass. The smallest sinks M < M© are drawn with twice the sink particle accretion radius 
race = 20 AU, the sinks with 1 Mq < M < 10 Mq with 5 race, and the sinks with M > 10 Mq with 10 race- 


next ~ 5 kyr, the bulk of sink formation occurs in a pair of tur¬ 
bulent, rotating, disky structures separated by ~ 2000 AU, as seen 
in the middle-left panel of Figure The structures subsequently 
merge into a larger such structure. Still later, a new pair of rotating 
disk structures is apparent. Additional sink formation occurs in two 
linear filamentary features extending in opposite directions, as seen 
in the middle panel of Figure Close encounters between sinks 
are common. Numerous sink ejections occur as a result of close en¬ 
counters and many-body interactions. The sink-gas system reaches 
approximate virial equilibrium relatively quickly, in ~ 5 kyr, after 
the onset of sink formation. 


Figure shows the evolution of the mass-weighted gas tem¬ 
perature and is otherwise identical to Figure Within ~ 8 kyr 
after the onset of sink formation, the dust temperature in the central 
~ 2000 AU begins to exceed ~ 100 K and the dust starts heating 
the gas. At the end of the simulation, the central ~ 3000 AU has an 
average gas temperature ~ 300 K, hot enough to significantly raise 
the thermal Jeans mass and inhibit new sink formation. 
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Figure 3. Same as Figurel^but for gas temperature. 


3.2 Sink Particle Formation and Growth 

We ran the simulation for ISkyr after the formation of the first 
sink. In this period, 30 sinks formed with a total mass of 81 Mq. 
For most of this time the median sink mass stayed constant at 
0.5 M 0 , though in the last ~ 3kyr, the median mass increased 
above 1 Mq . At the end of the simulation, the most massive sink 
was ^14 Mq and was undergoing Kelvin-Helmholtz contraction 
towards the main sequence. 

Figure [^displays a selection of sink properties as a function 
of time: the total accretion rate onto all sinks the total ac¬ 

cretion luminosity and intrinsic luminosity from all sinks, the total 
number of sinks, the total mass in sinks and the masses of individ¬ 
ual sink particles, the effective temperature of the four sinks that 


would grow to be the most massive ones by the end of the simula¬ 
tion, and the mean and median sink particle mass. 

Beginning 2 kyr after the formation of the first sink and 
until the end of the simulation, the total mass accretion rate onto all 
sinks grew from Mq yr“^ to almost Mq yr“^. In 

the final 2 kyr, there is indication of a slight downturn in 
though it is unclear if this is a genuine trend or a transient accretion 
rate fiuctuation. However, there is indeed a significant paucity of 
new sink formation in the final 4 kyr. 

In Figure we show the accretion rate as a function of 
time since birth for each of the 30 sinks that form in the sim¬ 
ulation. These accretion histories vary widely, but some general 
conclusions can be drawn. With very few exceptions and inde¬ 
pendent of the final sink mass, each sink accretes very slowly at 
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Figure 4. Sink particle properties as a function time, with the time coordinate shifted so that the first sink forms at f = 0. The top-left panel shows the total 
accretion rate onto sinks (blue) and the accretion rate onto the sink that grows the most massive by the end of the simulation (red). The top-right panel shows 
the total mass in sinks (red) and the masses of individual sinks (blue). The middle-left panel shows the total accretion luminosity (red) and total intrinsic 
luminosity (blue) for all sinks. It also shows, in green, what the total accretion luminosity would have been if the protostellar radii had been calculated using 
Equation instead of with MESA evolutionary tracks. The middle-right panel shows the effective temperatures of four sinks that are the most massive at the 
end of the simulation. The dashed line at Tgft = 10^ K is the approximate effective temperature limit above which a protostar would emit significant ionizing 
radiation. The bottom-left panel is the total number of sinks and the bottom-right panel shows the mean (blue) and median (red) sink mass. 


birth, M* < M© yr“^. The accretion rate rises to between 
10“^ M© yr“^ and 10“^ M© yr“^ in less than 1 kyr. We refer to 
this rise as the initial accretion rate peak. Such accretion rates are 
~ 2 — 20 times the nominal core accretion rate ^ cl/G expected 
in a marginally Jeans-unstable collapsing core with a temperature 
50K(e.g., |Shu|1977| ). 


The instantaneous accretion rate onto a given sink also ex¬ 
hibits strong temporal fluctuations, often by more than three orders 
of magnitude on time scales <C 1 kyr. Some of this variability is 
unphysical numerical discreteness of our operator-split scheme for 
accreting gas mass onto sinks. When M* is binned on 300 year time 
scales (red lines in Figure [^, short term variation is still present 
but with a lower amplitude. This variability in the accretion rate 
translates directly into variable accretion luminosity via Equation 
0. Observed embedded protostars typically have luminosities one 
to two orders of magnitude smaller than what would be naively 
expected from indirect estimates of accretion rates (e.g., |Kenyon| 
|et al.|p^^ [Evans et'd^|2009[ [Enoch et~n?]|2009| ). Highly vari¬ 
able, or episodic, accretion is a natural consequence of protostellar 
growth in a turbulent environment or from gravitationally unstable 
circumstellar disks. Variable accretion is one potential solution to 


the ‘protostellar luminosity problem’ (e.g., [Dunham & Vorobyov[ 

[ 201 ^ . 

Sinks that ultimately reach high mass, here arbitrarily de¬ 
fined as > 3M© (e.g., panels 1, 2, 3, 4, 5, 6, and 12 in 
Eigure [^, do so via extended periods of rapid accretion, with 
M* > M© yr“^ over > 10 kyr. Most of this accretion oc¬ 
curs in competition with other sinks and does not sample the mate¬ 
rial originally exclusively bound to, or even geometrically associ¬ 
ated with the sink. The lower panels of Eigureshow that a sink’s 
initial Jeans-unstable progenitor core is typically accreted entirely 
within ~ 3 kyr. At the end of the simulation, many massive sinks 
sustain high accretion rates, although we caution that our simula¬ 
tion does not include two principal mechanisms that could suppress 
accretion, protoionization and protostellar outflows. A modeling of 
protostellar outflows (e.g., [Cunningham et al.[201 l[[Eederrath et al. 


20141, photoionization (e.g., Peters et al.[[2010[ [Dale & Bonnell 


2011 1 , and a more accurate treatments of radiative transfer, would 


increase the physical realism, particularly with regards to the for¬ 
mation of massive stars and the driving of interstellar turbulence. 


Low mass sinks with M* < 1 M© (e.g., panels 8, 13, 14, 15, 
18, 22, 23, 26, and 27 in Eigure]^ typically exhibit sudden accre¬ 
tion rate termination, with the rate dropping below 10“^ M© yr“^ 
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Figure 5. Instantaneous accretion rates of individual sink particles as a function of time since the individual sink’s birth. The birth rank and the final mass at 
the end of the simulation Mf are shown in the top-left corner of each panel. The mass of the sink at the final simulation snapshot,, is also shown for each sink. 
The blue line shows the accretion rate extracted directly from the simulation and the red line is the accretion rate binned over a 300 year period. The vertical 
dashed green line in each panel shows the end of the simulation. Panels in which the accretion rate drops below 10“^ Mq yr“^ (e.g., panels 8, 13, 15, and 
23) show sinks that have experienced suppression or termination of accretion, usually after a close encounter with more massive sinks. 


and in some cases much lower. This can be attributed to the rel¬ 
atively compact star forming environment ( r\j 10'‘ AU) in which 
close encounters between sinks are commonplace. Many-body en¬ 
counters give large kinetic energy kicks to the lowest-mass bodies, 
placing them on fast or even unbound orbits. Such orbits sample 
lower density gas, which in turn reduces the Bondi-Hoyle accre¬ 
tion rate (Mbh oc pv~^). Many-body encounters and ejections of 
this type become progressively common as the crowding of the star 
forming complex increases. 

We observe an interesting correlation between birth rank and 
the initial peak accretion rate sinks achieve in their first ~ 1 kyr. 


The earliest sinks (panels 1 — 12) have typical initial peak accretion 
rates ~ 3 x 10“^ Mq yr“^ while later sinks (panels 17 — 30) have 
initial peak rates typically exceeding 10“^ M© yr“^. This can be 
understood as a consequence of stellar radiative feedback from ex¬ 
isting protostellar sources. Protostars heat the central star forming 
environment to an increasing degree and this raises the core accre¬ 
tion rate ~ cf/G oc in newly collapsing cores (see Figure]^. 
Many of these late-forming sinks with high initial peak accretion 
rates terminate accretion shortly after birth (e.g., panels 22, 23, 25, 
26, and 27 in Figure thanks to close encounters with massive 
sinks that are common in the dense star-forming complex. 
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We observe evidence for competitive accretion between the 
protostars in which accretion of initially unbound gas from a com¬ 
mon gas reservoir determines the final protostellar mass (e.g., |Bon^ 
|nell et al. 1 1998 |[200Ta| [Bate & Bonnell|2005| l. More massive stars 
not only have larger accretion cross sections, but they also segre¬ 
gate toward the cluster center, where the gas density, and thus accre¬ 
tion rate, are higher. Elementary analytical models of this ‘rich-get- 
richer’ scenario produce self-similar protostellar growth and pre¬ 
dict a scale-free, power-law mass distribution akin to the observed 
high-mass stellar IMF (e.g., |Bonnell et al.|2001b| l. Nevertheless, a 
predictive, realistic theory of the IMF taking into account competi¬ 
tive accretion remains lacking. Models not allowing for the dynam¬ 
ically complex competitive accretion, but rather assuming that stel¬ 
lar masses are imprinted by turbulent density fluctuations preceding 
star formation, have had more success in producing definitive pre¬ 
dictions for the IMF ( McKee & Tan|2002l |2003[ [Krumholz et al.| 


2005 [ IPadoan & Nordlunq|2002| |Hennebelle & Chabrier||2008l 


Hopkins|[20I2|l, in spite of the inevitable role of competitive ac¬ 


cretion. Our sinks continue growing well after accreting the initial 
gravitationally bound ‘core’ (Figureby accreting in a relatively 
crowded star-forming environment (Figurej^. The sinks most mas¬ 
sive at the end of the simulation acquire the bulk of their mass in the 
latter regime. This underscores the importance of the competitive- 
accretion-like mode of gas accretion, at least in the specific regime 
explored here. 

We should note that the initial conditions play a crucial role 
for the mode of gas accretion and stellar growth, as has been em¬ 
phasized by numerous authors (e.g., [Krumholz et al.|2005[[Bonnell| 
|& Bate|2006[|Girichidis et al.|2011| ). In particular, simulations ini¬ 
tialized with decaying turbulence tend to produce a coherent, cen¬ 
tral collapse, ideal conditions for competitive accretion. Those ini¬ 
tialized from driven turbulence, in which small scale turbulence is 
continually replenished from larger scales, seem instead to fix stel¬ 
lar masses through fragmentation of turbulent structures into self- 
gravitating cores. The initial conditions of the present simulation, 
which were extracted from a parent cosmological simulation at a 
site of thermal instability, can be regarded as belonging in the for¬ 
mer category, containing turbulence initially excited by thermal in¬ 
stability and then decaying, or growing if amplified by gravitational 
compression (e.g., [Robertson & Goldreich|2012j |. 


3.3 Sink Particle Mass Function 

In Figure we show the mass function of sink particles between 
0.1 Mq and 20 M© in three representative simulation snapshots. 
The mass function evolves considerably over ~ 10 kyr, particu¬ 
larly in view of the late formation of > 3 Mq sinks. In blue, we 
also show the mass function of sinks that have ‘stopped grow¬ 
ing’, here defined as the sinks with accretion rates remaining be¬ 
low 10~^ Mq yr“^ for the final 0.5 kyr of the simulation. These 
are typically the sinks that have had a strong dynamical encounter 
with higher mass sinks. All but two of the non-accreting sinks have 
masses below 1 Mq at the end of the simulation. 

Above ~ 3Mq, the sink mass function develops hints of a 
power-law-like tail toward higher masses. Since there are only 12 
sinks with such masses, the statistics is insufficient to perform esti¬ 
mation of the power-law slope. We also caution that the sink mass 
function has not converged at high masses. Had we run the sim¬ 
ulation longer, the high-mass tail would have undoubtedly grown, 
tempered, perhaps, by the effects of photoionization that would be 
essential to include in the high-mass regime. 

In the early-results version of the present simulation 



Figure 6. Sink particle mass function at 7.6 kyr (top panel), 11 kyr (middle 
panel), and 18 kyr (bottom panel) after the onset of sink formation. The red 
histograms are for all sink particles and the blue histograms are only for 
non-accreting sinks, namely those with instantaneous accretion rates below 
10“^ Mq yr“^ in the last 500 yr of the simulation. The black and green 
straight li nes indicate pawer -low slopes of o = 2.35 jSalpeter|1955) and 
a = 1.3 jCeha et al.|2013} where dN/dM^ oc . The solid portion 
of the Geha et al. line, between 0.5 Mq and 0.8 Mq, indicates the narrow 
stellar mass range in which observations constrain the IMF in two ultra-faint 
dwarf spheroidal satellites. 


( [Safranek-Shrader et al.[[^14a[ ), which had not been run long 
enough to form high mass > 3 M© sinks, we tentatively at¬ 
tempted to constrain the power-law slope of the sink mass func¬ 
tion in the intermediate mass range 0.1 M© < M* < 3M©. 
We obtained a maximum likelihood slope of a 1.3 such that 
dN/dM^ oc M~^. Perhaps coincidentally, this slope was consis¬ 
tent with the stellar mass function slope a = 1.3 measured in the 
narrow stellar mass range 0.5 M© < M* <0.8 M© in ultra-faint 
dwarf (UFD) satellite galaxies Leo IV and Hercules ( [Geha et 
[2013[ ). In the present, extended simulation, the mass function in the 
range 0.1 M© < M* < 3M© seems fiat with a ^ 1, similar 
to, and perhaps even shallower than reported in [Geha et al.[ ( [2013] ) 
and [Safranek-Shrader et al.[ ( |2014a[ ), albeit with bin-to-bin statisti¬ 
cal fiuctuations, including an apparent, but statistically insignificant 
trough at M* ~ 2 M©. 

We reiterate, however, that all sinks above 1 M© are still ac¬ 
creting at the end of the simulation (Fig.[^ and the sink mass func¬ 
tion is far from convergence. The results reported here should be 
construed as diagnostic of a trend in the protostellar mass function 
assembly process, rather than predictive of the final stellar IMF. It 
is also important to emphasize that the sink particle mass need not 
directly translate in the final stellar mass because only a fraction of 
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the mass accreted by a sink particle is incorporated into the final 
star. For example, mass loss to protostellar outfiows is expected re¬ 
duce the core-to-star conversion efficiency by ^ 50% (|Matzner &| 
|McKee|2000||Alves et al.|2007l[Enoch et aL|2008| t. 

There are other effects that could suppress gas accretion and 
inhibit further evolution of the mass function. HII regions produced 
by massive stars disperse natal molecular clouds and shut off accre¬ 
tion (e.g., |Whitwor4L|1979[[Matzner|20'0^|Krumholz et al.|2006| l. 
This could halt the evolution of the stellar mass spectrum, perhaps 
not long after the end of the present simulation. Indeed, as shown 
in the middle-right panel of Figure]^ the most massive sink has an 
effective temperature of Tgff 34, 000 K, hot enough to produce 
ample ionizing luminosity. However, an intense, sustained gas in- 
fiow can curtail an HII region ( |Walmsley|1995[[Keto|2003| l. Also, 
filamentary and disky structures that form in turbulent star-forming 
environments tend to resist photoionization (e.g., [Peters et al.|2010[ 
[Dale & Bonnell|2011| |. Massive stars preferentially ionize low den¬ 
sity gas, allowing the bulk of accretion to continue through low- 
filling-factor channels. 

3.4 Protostellar Evolution and the Heating of Dust and Gas 

For the first ^12 kyr of sink growth, accretion luminosity dom¬ 
inated protostellar radiative output. Intrinsic luminosity took over 
when the first protostar, at 7 M©, underwent core contraction, 
transitioning from the convective Hayashi track to the radiative 
Henyey track. Its radius, effective temperature, and intrinsic lumi¬ 
nosity increased significantly and abruptly, as is clear in Figure 
and Equation This intrinsic brightening is also evident in the 
middle-right panel of Figure where the effective temperatures 
of the most massive sinks rapidly increased from 5000 K to 
> 25, 000 K in a time span of only ~ 1 kyr. At the end of the 
simulation, the most massive sink has an effective temperature of 
3.3 X 10^ K, hot enough to emit a significant ionizing luminos¬ 
ity, an effect we do not include in this study. The total protostellar 
luminosity is nearly equally divided between the accretion lumi¬ 
nosity from all sinks combined and the intrinsic luminosity from 
the three sinks that have begun hydrogen burning. 

Thermodynamics modulates fragmentation in the course of 
gas collapse. In Figure we explore the thermodynamic evo¬ 
lution in the density-temperature phase plane. We plot gas and 
dust temperature against density in four representative simulation 
snapshots. At the beginning of the excised simulation at redshift 
z 13.8, the initial gas temperature equals the CMB temperature 
Tcmb = 40 K, while the peak gas density is nn ~ 10^ cm“^. As 
gas collapse proceeds, two effects act to gradually raise the temper¬ 
ature. First, [CII] and [O I] cooling becomes less effective as these 
lines become optically thick, both intrinsically and to the dust con¬ 
tinuum. Second, the gas is slightly heated by the exothermic for¬ 
mation of H 2 on the surfaces of dust grains and to a lesser extent in 
three-body reactions. This is evident in Figure[7^, showing the state 
of the gas after 9.6 kyr. However, the gas temperature drops back 
to Tcmb soon thereafter, when at density nn ~ 10® cm“®, 
gas becomes thermally coupled to dust. In Figure |7]3, gas be¬ 
comes optically thick to the continuum at nu ~ 10”cm“® and 
begins to evolve adiabatically until sink formation is possible at 
nn = 4 X 10^^ cm“®. Protostellar radiative heating of dust grains 
alters this basic picture somewhat, by increasing the thermal Jeans 
mass in proximate collapsing cores. This can be seen in Figure]^ 
recorded after the first sink has grown to ~ 0.5 M©. Finally, Fig- 
ure[^ shows the thermodynamic state at the end of the simulation, 
18 kyr after the first sink formed. The heating of gas via radiative 



102 ^q6 ^q8 ^qIO ^q12 

Hh [cm-^] 



102 10^ 10® 10® 10’° 10'^ 




Hh [cm-^] 

Figure 7. Thermodynamic evolution of the gas and dust. Colored cells rep¬ 
resent the amount of gas in density-temperature cells, with red representing 
the highest gas mass per cell. The dust temperature is overplotted with un¬ 
filled rectangles. Dashed lines indicate representative values of the Jeans 
mass. The panels show the state: (a) 9.6 kyr after the start of the simulation 
when the gas begins to depart from isothermal collapse, (b) at the instance 
of the first sink particle formation, (c) 1.5 kyr later, when protostellar ra¬ 
diative heating of dust starts to refiect in the gas temperature, and (d) at the 
end of the simulation. 
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heating of dust is strong at densities nu > 10^ cm“^ at which the 
gas temperature increases with increasing density. 

The early thermodynamic evolution of the first protostellar 
core undergoing gravitational collapse in Fig ure [^,b fits the pre- 
dictions of idealized, one-z one models (e.g.,|Oniukai et aL|2005[ 

[Schneider & Omukai|2010[[bmukai et al.|2010| |. By following the 

gravitational collapse and thermodynamic evolution of gas in one 
spatial dimension, jOmukai et aL| ( [20T0] ) estimated that accretion lu¬ 
minosity does not suppress fragmentation for gas metallicities be¬ 
low 10“^ Zq, but only considered the impact on the birth core. 
Figure [^,d, however, shows that as protostars grow, radiative heat¬ 
ing of dust grains produces significant departures from one-zone 
collapse models (also see jDopcke et al.|2013j t. Late accretion onto 
sinks and the formation of new sinks in sites pre-heated by pro¬ 
tostellar radiation are significantly modified by dust heating. This 
highlights the important role that protostellar radiative feedback 
plays in modulating the star formation process, even at sub-solar 
metallicities. 


10 ^ Zq at redshift z = 13.8, one-zone models predict a charac¬ 
teristic fragmentation mass of ~ 1 M© (e.g., jOmukai et al. 120051 
[Schneider & Omukai|2010[[Schneider et al.|2012| l. This mass stems 

from the onset of dust-gas coupling and departure from isothermal- 
ity that should set in at densities nu ~ 10® cm“®, in good agree¬ 
ment with our findings. It supports the hypothesis that dust grains, 
and not metal fine-structure lines, are responsible for moderating 
the Pop III to Pop II star formation transition. Metal line cooling, 
however, plays a dominant role in regulatin g fragmentation at more 
moderate densities nu ~ 10® — 10® cm“® ( Safranek-Shrader et al. 
[2014b[ ). 


If indeed the characteristic stellar mass variation can be mod¬ 
eled as in [Schneider & OmuHi] ( [2010[ ), the relationship between 
characteristic mass, redshift, metallicity, and dust content can be 
used as a probe of ancient stellar populations via dwarf galaxy ar¬ 
chaeology (e.g., [Frebel & Bromm[|2012[ l. Further hydrodynamic 
studies surveying a wider parameter space are necessary to validate 
the applicability of one-zone models to quantifying star formation 
in realistic, multi-dimensional, turbulent environments. 


3.5 The Imprinting of a Characteristic Mass 

The thermodynamic behavior discussed in the previous sub-section 
seems to plays a central role in imprinting a characteristic mass on 
the fragmentation process. Observations have shown that the char¬ 
acteristic stellar mass scale, i.e., the median stellar mass, is remark¬ 
ably constant in the Galaxy across a wide range of stellar environ¬ 
ments (e.g., [Bastian et al.|2010[ |. A median mass that varies with 
metallicity, redshift, halo mass, or other factors would have im¬ 
plications for any remotely detailed modeling of reionization and 
chemical enrichment and especially for the interpretation of extra- 
galactic observations. In the present simulation, there is a nearly 
time-invariant median sink particle mass of M^c ~ 0.5 M©. Sinks 
with higher masses seem to grow competitively, by accreting from 
a gas reservoir shared with other sinks. On the other hand, below 
the median mass, momentum kicks acquired in close and many- 
body encounters terminate accretion, leaving the sink with its ini¬ 
tial bound core mass (see Fig. and [^. What is the physics that 
selects M* as the characteristic mass scale? 

Since the equations of ideal, isothermal, self-gravitating 
(magneto-) hydrodynamics are scale-free (e.g., [Krumholz|[2014[ |, 
explanations for a characteristic mass scale require a departure 
from isothermality]^ This departure is generally attributed to ei¬ 
ther protostellar feedback processes, such as protostellar radiative 
heating of gas that couples thermodynamics to star formation ( [Bate[ 
[2009[ [Krumholz[[201 1[ ), or chemical, radiative, etc., processes in¬ 
trinsic to the gas that select specific temperature and density scales 
(e.g., [Larson|2003[ ). In the remainder of this sub-section, we focus 
on the latter. Idealized, one-zone models for the gravitational col¬ 
lapse and fragmentation of metal-poor clouds have identified the 
redshift (which sets the CMB temperature fioor), gas metallicity, 
and dust-to-gas ratio as the parameters with the strongest influence 
on the characteristic fragment ation mass (e.g.,[Omukai et al.|2005[ 
[Safranek-Shrader et al.|2010[ [Schneider & Omukai|2010| ). These 

studies, however, did not explore many other effects such as radia¬ 
tive feedback, magnetic fields, or supersonic turbulence. 

For the system in our simulation with metallicity Z = 


^ Alternatively, a characteristic mass scale could also be imprinted by a 
fine tuning of initial conditions by a process preceding, and external to star 
formation. We do not find this hypothesis attractive. 


4 DISCUSSION AND CONCLUSIONS 


[Safranek-Shrader et al.[ p014b[ ) and the present companion paper 
taken together trace the process of low-metallicity, primordial star 
formation from the initial conditions imprinted in cosmic density 
fiuctuations to the for mation of individual protostars. In [SafraneE^ 


Shrader et al. 


p014bt we simulated a 1 comoving Mpc® cosmolog¬ 


ical volume beginning at z = 145 until a halo with Tvir = 10 K 
has virialized at z = 13.8. At that point, the initially metal-free 
gas was assigned a non-zero metallicity of 10' ■'^0, crudely mod¬ 
eling chemical enrichment by preceding Pop III stars. The dens¬ 
est, metal-enriched gas cooled isobarically via [CII] and [O I] line 
emission to Tomb = 40 K. The cold gas fragmented into sev¬ 
eral few-hundred-solar-mass clumps. The present simulation was 
initialized by excising one of these fragmentary clumps, contin¬ 
uing the simulation to densities nu > 10^^ cm“®, and inserting 
sink particles to track individual protostars in the clump. Over the 
course of 18 kyr of star formation, 30 sink particles formed with 
a combined mass of 81 M©. The individual sink masses ranged 
from 0.05 M© to 14.4 M©. Massive stars had a strong radiative 
and dynamical impact on the star formation process, stunting the 
growth of lower mass protostars and suppressing the formation of 
new ones. 

The present simulation overlaps with and extends the early run 
reported in [Safranek-Shrader et al.[ ( [2014a[ ), hence it merits exam¬ 
ining any potential differences between the two runs. After 7 kyr 
of protostellar formation, the HEAT simulation of [ Safranek- Shrader[ 
[et al.[ ( [2014a| ) formed 37 sinks with a total mass of 15M©. In 
the same time period, the present simulation formed only 12 sinks 
with a total mass of 15 M©. While the star formation efficiencies 
match, a factor of 3 fewer sinks in the present simulation de¬ 
serves scrutiny. Two differences between these otherwise identical 
simulations contribute to this discrepancy. The present simulation 
was run at a slightly lower resolution, with a sink particle accretion 
radius of race = 20 AU compared to 10 AU in [Safranek-Shrader[ 
[et al.[ ( [2014a[ t. The new reduced resolution is still sufficient to al¬ 
low the gas to cross the opacity limit for fragmentation (Equation 
1^ preceding sink particle creation. A more probable cause for the 
discrepancy is different subgrid prescriptions for protostellar evolu¬ 
tion, as disc ussed in Section [2.4[ Com pared with the approximate 
treatment in [Safranek-Shrader et al.[ ( [2014a[ ) that assumed quasi- 
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spherical, photosphere-blanketing accretion, the MESA generated 
protostellar evolutionary models assume thin disk accretion and 
yield much smaller protostelar radii, and thus much larger accre¬ 
tion luminosities. Indeed, the green curve in the middle-left panel 
of Figured shows that if the expression from | Stabler et al.| ( |1986| ) 
(our Eqn. were used instead of the MESA tracks, the total ac¬ 
cretion luminosity would be roughly an order of magnitude lower, 
except immediately after sink insertion. 

The onset of metal line cooling in the parent virialized cos¬ 
mological object triggered thermal instability that allowed gas to 
cool quasi-isobarically to the CMB temperature. At this point, frag¬ 
mentation set in, so in |Safranek-Shrader et al.| ( |^14b| ), coarse sink 
particles were utilized to follow the evolution of the fragments for 
~ 4 Myr. The coarse sinks did not represent individual protostars 
but instead pre-stellar clumps (e.g., |Bergin & Tafalla|2007l l, each 
likely to fragment into multiple protostars. In the cosmological run 
with a metallicity matching the one in the present simulation, 11 
such clumps formed separated by 100 — 500 kyr intervals. We re¬ 
iterate that the present paper investigates the formation of stars in 
one of these clumps, specifically the first one to form, and follows 
its internal evolution for 60 kyr. 

If our results are representative of how stars form in all clumps 
and we can calibrate their star formation efficiencies to that of the 
first clump, we estimate that our canonical metal-enriched atomic 
cooling halo forms a stellar mass of ~ 300 M© in all clumps after 
~ 10 kyr of stellar growth in each clump. If the entire mass of 
each pre-stellar clump is converted into stars over a longer period 
of time, the final stellar mass is ~ 3000 M©. Even the latter, more 
optimistic estimate of the star formation yield implies a relatively 
low halo-wide star formation efficiency of M*,tot/Afb ~ 6x10“^, 
where is the total stellar mass and Mb is the total baryonic 

mass in the halo. 

These calculations allow us to assess the prospect that star 
clusters of this kind can be detected with next-generation tele¬ 
scopes, such as the JWST. Eollowing |Pawlik et al. (|2011|), we 
first use the stellar population synthesis models of |Schaerer||2003| ) 
to estimate the initial luminosity L(Ha) in the Ha line and the 
luminosity L(Hel640) in the Hell line, as well as recombina¬ 
tion line and the initial ultraviolet (UV) continuum luminosity per 
unit frequency L(UV1500) at the restframe wavelength 1500 A 
as produced by an instantaneous, metallicity Z = 10“^ Zq star- 
burst forming a 3000 M© star cluster. This yields L(Ha) = 3 x 
10^®ergs-\L(Hel640) = 6x10^^ ergs-\andL(UV1500) = 

2 X 10^^ ergs“^ Hz“^. Eor a spatially unresolved source, we 

can translate these luminosities into observed fiux densities us¬ 
ing, e.g.. Equations (4) and (5) of |Pawlik et al.| ( [^ll| ). Eor the 
combined stellar and nebular UV continuum intensities, we find 
/iy(UV1500) 10“^ nJy. Similarly, assuming a spectral resolu¬ 

tion of R = A/A A = 3000 for Ha and R = 1000 for Hel640, we 
find fu{Ha) ^ 1 nJy and /iy(Hel640) ^ 2 x 10“^ nJy. Given a 
source redshift of z = 13.8, JWST will detect the redshifted UV 
continuum at Ao = 2.2 /xm with the NIRCam instrument, Ha at 
Ao = 9.8 /xm with the MIRI spectrograph, and the He II recombi¬ 
nation line at Ao = 2.4 /xm with the NIRSpec spectrograph. With 
a 10^ s exposure, the sensitivities of these three instruments are 10, 

3 X 10^, and 9 x 10^ nJy, respectively, to achieve a signal-to-noise 
ratio of 10, clearly many orders of magnitude higher than the antic¬ 
ipated fiuxes. 

Even assuming (incorrectly) that the IME would grow to be 
top-heavy, which would substantially increase the stellar and neb¬ 
ular luminosities—particularly in the He1640 line—it is clear that 
the cluster does not fall within the deep, broadband or spectroscopic 


detection limits of JWST (also discussed in, e.g., [Gardner et al.] 
|2006| [Johnson et al.|2009[[Pawlik et al.|2011[ [Smith et al.|2014| |. 

The cluster may be magnified by foreground gravitational lenses 
(e.g., [Zackrisson et al.[2012[[2014[ ), but unfortunately, even extreme 
magnifications (/x 100) are insufficient for a direct JWST detec¬ 

tion. The high-redshift, metal-enriched stellar systems that will be 
detected with JWST will not be the very first metal-enriched stel¬ 
lar systems, but more evolved galaxies that have almost certainly 
been polluted by multiple generations of supernovae. These more 
evolved systems likely reside in much h 
with virial masses Mvir ^ 10^ M© (e.g.. 

These considerations indicate that as 
III and the earliest Pop II star formation must rely on alternatives 
to direct point source detection at a high redshift. A promising ap¬ 
proach involves chemical mapping of metal-poor stars in the Galac¬ 
tic halo and UED satellite galaxies. Owing to their very low metal- 
licities, chemical heterogeneity, and prevalence of old stars, these 
systems are considered to be relics of the first galaxies ([Ricotti[ 

& Gnedin||2005||Bovill & Ricotti||2009l IFrebel & Bromm||2012| 

Brown et al.||2013| |2014) . The stellar system we simulated here 
shares physical characteristics with the faintest UEDs—Bootes II, 
Segue I and II, and Willman I—in terms the dark matter mass col¬ 
located with the stars, the estimated stellar mass, and metallicity. 
Simulations tracking the first metals to their progenitor Pop III 
supernovae, following the subsequent chemical transport at high 
numerical resolution, and allowing for fine-grained chemical het¬ 
erogeneity, are necessary to further elucidate the relation of the 
first metal-enriched star forming systems to UEDs (see [Ritter et al.| 
\20l7\\2m^ . 


jirger dark matter h aloes. 


Pawlik et al.|201l| ) 


tronomical probes of Pop 


Published simulations of the onset of metal-enriched star for¬ 
mation in the first galaxies ( Wise & Abel|2008[[Greif et al.|2010 


Maio et al.|20I0|[Aykutalp & Spaansl201I [Wise et al.|20I2[ Mu-[ 

ratov et al.|2013[ [Johnson et al.|20I3|Veon et al.|2014[[201^ have 


yet to resolve the differentiation of gas clumps into individual pro¬ 
tostars and could not derive the star formation efficiency or the 
shape of the IMP from first principles. Our simulation was specifi¬ 
cally designed to return predictions of these key variables. While a 
single simulation does not allow us to explore statistical variation, 
it demonstrates that computationally expensive “direct”, yet fully 
cosmological, simulations of star formation can be used to cali¬ 
brate subgrid models that can then be used in coarser “large-eddy” 
simulations of star formation in much larger cosmological volumes 
and across longer redshift intervals. This is particularly crucial for 
modeling the early stages of cosmic reionization, currently plagued 
by substantial uncertainties regarding the star formation efficiency 
in the smallest galaxies. 
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